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ABSTRACT 

Upcoming facilities such as the Herschel Space Observatory or ALMA will deliver a wealth of molec- 
ular line observations of young stellar objects (YSOs). Based on line fluxes, chemical abundances can 
then be estimated by radiative transfer calculations. To derive physical properties from abundances, 
the chemical network needs to be modeled and fitted to the observations. This modeling process is 
however computationally exceedingly demanding, particularly if in addition to density and tempera- 
ture, far UV (FUV) irradiation. X-rays, and multi-dimensional geometry have to be considered. 
We develop a fast tool, suitable for various applications of chemical modeling in YSOs. A grid of 
the chemical composition of the gas having a density, temperature, FUV irradiation and X-ray flux 
is pre-calculated as a function of time. A specific interpolation approach is developed to reduce the 
database to a feasible size. Published models of AFGL 2591 arc used to verify the accuracy of the 
method. A second benchmark test is carried out for FUV sensitive molecules. 

The novel method for chemical modeling is more than 250,000 times faster than direct modeling and 
agrees within a mean factor of 1.35. The tool is distributed for public use. Main applications are (i) 
fitting physical parameters to observed molecular line fluxes and (ii) derive chemical abundances for 
2D and 3D models. They will be presented in two future publications of this series. 
In the course of devloping the method, the chemical evolution is explored: We find that X-ray chem- 
istry in envelopes of YSOs can be reproduced by means of an enhanced cosmic-ray ionization rate 
with deviations less than 25%, having the observational consequence that molecular tracers for X- 
rays are hard to distinguish from cosmic ray ionization tracers. We provide the detailed prescription 
to implement this total ionization rate approach in any chemical model. We further find that the 
abundance of CH+ in low-density gas with high ionization can be enhanced by the recombination of 
doubly ionized carbon (C"'"'") and suggest a new value for the initial abundance of the main sulphur 
carrier in the hot-core. 

Subject headings: molecular processes - methods: numerical - astronomical data bases: miscellaneous 
- stars: formation - ISM: molecules - methods: data analysis 



1. INTRODUCTION 

There is an interesting but little explored phase in 
star formation when the cloud core collapses, but the 
protostar is still deeply embedded. In this phase, the 
temperature of the interior envelope exceeds 100 K, 
outflows are observed, protoplanetary disks form, and 
the protostar begins to radiate in UV and X-rays. 
These physical processes can best be observed in low 
temperature lines, some atomic but most molecular. 
A wealth of new line data of young stellar objects 
(YSOs) will be available soon by upcoming facilities 
like the Herschel Space Observatory or the Atacama 
Large Millimeter Array (ALMA, e.g. van Dishoeck & 
J0rgensen 2008 for a recent review). Two steps are 
necessary to constrain physical parameters by molecular 
line observations: (i) Radiative transfer calculations 
are applied to estimate the abundance of the observed 
molecule, and (ii) modeling the chemical network 
relates the abundance to physical parameters such as 
age, density, temperature, far UV (FUV) and X-ray 
irradiation by the protostar. Although the present 
knowledge of these networks is still limited by unknown 
reactions on grain surfaces, in some cases it is possible 
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to derive chemical abundances from a set of initial 
physical parameters. The abundances form the bases for 
radiative transfer calculations. The resulting line fluxes 
can then be compared to the observed fluxes. Finally, 
the input physical parameters for the chemical modeling 
are changed until the derived line fluxes fit the observed 
ones, for instance producing a minimum in a test. 
Chemical modeling simulates a network of molecular 
species reacting with each other (e.g. Doty et al. 2002). 
Contrary to an atomic gas, molecular abundances are 
not conserved, but depend on local physical parameters. 
In chemical simulations, the abundances of different 
species are related to the change of one species in a 
set of non-linear coupled differential equations. Solving 
these equations is a time-consuming task for a chemical 
network consisting of several hundred species connected 
by thousands of reactions. 

While early astrochemical models of envelopes of 
YSOs assumed a fixed density and temperature (e.g. 
Leung et al. 1984), more recent models considered 
space-dependent chemistry (e.g. Caselli et al. 1993 
or Millar et al. 1997b). So far, chemical models have 
only been fitted to observations under the assumption 
of spherically symmetric, one-dimensional physical 
parameters (Doty et al. 2002; Stauber et al. 2005). This 
is a questionable assumption as YSOs have inherently 
2D or 3D geometries with outflow cavities, disks etc. For 
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such geometries the number of ceUs for which chemical 
modehng has to be performed increases from about one 
hundred to 10^ and more. Fitting requires repetition 
of many similar calculations. It becomes inefficient and 
beyond present resources for 2D or 3D geometries and 
for a large sample of sources. 

In a dynamical situation, the "chemical age" (evolu- 
tion time since the start from an initial composition) 
may be different from the age of the YSO. The chemistry 
in a gravitationally collapsing envelope has been studied 
by Rodgers & Charnley (2003) and Lee et al. (2004). 
Doty et al. (2006) have constructed a model including 
infall as well as the evolution of the central stellar 
object. They conclude that the chemical evolution of 
an ice-evaporated hot core appears to be younger than 
the time since the formation of the YSO by a factor 
of 4 to 10. "Pseudo time-dependent" models keep the 
physical structure constant with time, while the chem- 
ical abundances evolve. Such models nevertheless are 
reasonable approximations if the chemical composition 
evolves rapidly compared to the change of the relevant 
physical parameters, like for example in regions with 
strong FUV irradiation. 

Undoubtedly, the ionization of molecular hydrogen, 
H2, followed by fast ion- molecule reactions involves ex- 
tremely rapid processes. Ionized molecules are a driving 
force of the chemistry in the envelope of YSOs (Dalgarno 
2006). Cosmic rays account for ionizations mainly in the 
outer part of the envelope. In the inner part, high-energy 
radiation from the YSO may be the dominant source 
of ionization. Accretion phenomena and shocks may 
account for FUV radiation in low-mass stars, while 
hot high-mass stars emit copious FUV photons mostly 
from their hot photosphere. Direct photoionization 
through FUV photons leads to a particular chemistry 
as observed in photo dissociation regions (PDRs). FUV 
radiation is attenuated by dust (e.g. Montmerle 2001). 
Unless a low-density region allows FUV radiation to 
escape from the innermost region of the envelope, its 
influence on the chemistry is restricted to a small volume 
behind the irradiated surface (Stauber et al. 2004). The 
observations of Stauber et al. (2007) revealed a much 
larger amount of several FUV sensitive molecules than 
predicted by their spherically symmetric models. These 
authors suggested outflows acting as paths for FUV 
photons to escape and penetrate the high-density gas in 
the border region between the outflow and the envelope. 
A multidimensional geometry is needed to model these 
"outflow- walls" . 

In addition to FUV, young stars are known to be 
strong emitters of X-rays with luminosities up to 
10^^ erg s~^ in the 1 - 100 keV band (e.g. Preibisch & 
Feigelson 2005). Magnetic activity is believed to be the 
origin of the X-rays, but the exact mechanism remains 
unclear. X-ray photons have a smaller absorption 
cross-section than FUV (oc A"^) and penetrate deeper 
into dense material than FUV. X-ray ionization thus 
decreases mostly by geometric dilution, cx with 
distance r to the source. Various studies on the influence 
of X-ray irradiation on the chemistry of molecular clouds 
have been carried out (Krolik & Kallman 1983, Maloney 



et al. 1996, Lepp & Dalgarno 1996, Yan 1997, Stauber 
et al. 2005). For low- mass YSOs, the presence of X-rays 
has implications on the disc and thus planet formation 
(e.g. Ilgner & Nelson 2008; Meijerink et al. 2008). In an 
early stage of star formation, FUV and X-ray radiation 
are absorbed by the large opacity of the envelope and are 
therefore not directly observable. The point of evolution 
when FUV and X-ray activity sets in is unknown. 

The goal of our work is to provide a fast and simple 
method for chemical modeling. A large set of pre- 
calculated abundances for different ages and physical 
conditions (density, temperature, etc.) is used to quickly 
interpolate the chemical abundances for individual 
conditions in a model with physical conditions that 
change with position. The main problem is to reduce the 
high-dimensional space of physical parameter to a size 
that fits current hardware resources, but keep the accu- 
racy of the interpolation at an acceptable level. Based 
on this new method for chemical modeling, observed 
data can be quickly fitted to physical parameters such as 
the chemical age or the X-ray luminosity. Also detailed 
2D or 3D models may be constructed quickly allowing 
to study the influence of the geometry on observable 
parameters (Bruderer et al., in prep.) 

This paper is organized as follows: In the first section, 
we describe the grid approach. We discuss the relevant 
parameters for the chemical composition and the simi- 
larity between X-ray and cosmic ray induced chemistry. 
In Section 3, the chemical model (network and reaction 
rates) used for this work is briefly described. Two 
benchmark tests on realistic applications are carried 
out and discussed in the following section. Section 5 
describes the application of the method. In the following 
parts of this series of papers, we will present a detailed 
multidimensional chemical model of a high-mass star 
forming region applied to the enigmatic C0+ molecule 
and a parameter study of 2D effects on the interpretation 
of line observations as expected to be observed with 
Herschel and ALMA. 

2. A GRID OF CHEMICAL MODELS 

Chemical models solve for the molecular abundances 
using the rate equations 

^ = E • "(^■) + E Kj,k ■ n{j) ■ n{k) + 5(z) (1) 

which relate the temporal evolution of a species, labeled 
by i, to the number density n{j), n{k) [cm~'^] of some 
species j and k. The constants of proportionality fc^j 
[s^^] and k[ ^ j, [cm^ s~"^] depend on physical properties 
like temperature, cosmic-ray ionization rate or FUV 
flux. Some applications require an additional source 
term S{i) to account e.g. for molecules evaporating from 
dust grains or spatial flows. 

Many of the reactions relevant in the interstellar envi- 
ronment at low temperature are difficult or impossible to 
measure in terrestrial conditions. In standard networks 
for astrochemical applications, a majority of reaction 
rates are thus only known to within a factor of two 
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(Woodall et al. 2007). Furthermore, uncertainties in the 
rate cocfBcicnts may grow due to the non-linear nature 
of the rate equations (Eq. 1). Studies on the quantita- 
tive uncertainty of the chemical abundances have been 
carried out only recently by Wakelam et al. (2005, 2006) 
and Vasyunin et al. (2008). Considering their results, it 
seems reasonable to require the interpolated abundances 
to agree within a factor of two with the fully calculated 
values. 

The basis of our method is to interpolate abundances 
from a grid of physical parameters relevant for the 
chemical composition. These physical parameters 
form a multidimensional space, where the dimensions 
corresponds to temperature, density, cosmic-ray ion- 
ization rate and several parameters for the impact of 
high-energetic radiation (X-rays and FUV irradiation). 
For each grid point and all molecules in the chemical 
network, the temporal evolution is stored in a database. 
Therefore, time enters as an additional dimension. In 
order to explore a wide range in parameter space, we 
assume a rectangular grid: For every grid point along 
one axis, all possible combinations of grid points along 
all other axes are calculated. The total number of 
models is obtained by multiplying the number of grid 
points along each dimension. Thus, it is crucial to keep 
the number of dimensions and grid points as low as 
possible, while it must still be sufficient to guarantee the 
required interpolation accuracy. 

As an example for problems to be met by the inter- 
polation, we show the water fractional abundance with 
respect to H2 in Fig. la: A parcel of gas and dust at 
a gas density between 10** cm~^ and 10* cm~^ ir- 
radiated by X-rays has been modeled in a similar way 
as Stauber et al. (2007). For this figure, the fractional 
abundances are given at the time (the so called chemical 
age) of 5 X 10'* yr according to the age of the high-mass 
protostellar object AFGL 2591 (Stauber et al. 2005). 

The presentation of Fig. 1 with logarithmic scales in 

both dimensions suggests to use a linear interpolation of 
the fractional abundance in log-log space. The straight 
(red/grey) linos in the figure show interpolated fractional 
abundances. The temperatures marked by arrows at 
the top of the figure act as interpolation points. Two 
temperatures close to 100 K have been used, since water 
is assumed to evaporate at this temperature from dust 
grains. For an accurate interpolation in the temperature 
range between f»200 and 400 K, a couple of interpolation 
points is required. The reaction OH -|- H2 — > H2O -|- 
H with a high activation energy becomes important in 
this temperature range, leading to little net variation 
in the water abundance at T > 400 K, while at lower 
temperature H2O is destroyed by X-rays (Stauber et al. 
2006). 

2.1. TemperatAire and density dependence 

The dependence of abundance on parameters can be 
much simpler than the temperature sensitivity of water: 
The fractional abundances of H3, HCO"'", H2O and 
H2CO depending on the total gas density are shown in 
Fig. lb. For that plot, a temperature slightly above 



the water evaporation temperature of 100 K has been 
assumed. No X-ray or FUV irradiation is considered, 
but cosmic rays ionize H2 and account for the production 
of Hg". A cosmic-ray ionization rate Cc.r. = 4 x 10~^^ 
s-i, found by Doty et al. (2004) in IRAS 16293-2422 is 
used for this plot. Higher values than the "standard" 
value of a few times 10^^^ s^^ have also been suggested 
by van der Tak & van Dishoeck (2000). This shows the 
necessity to include the cosmic-ray ionization rate as a 
dimension for the interpolation approach. The rate of 
ionizations by cosmic ray particles and thus the density 
of Hg is independent of the H2 density. Therefore, the 
fractional abundance of H3 is approximately inversely 
proportional to the density. The same effect is found 
for HCO"*", produced by the reaction Hj^ + CO 
HCO"*" -I- H2. In contrast, water and formaldehyde are 
not directly related to the ionization rate and their 
fractional abundances does not vary more than an order 
of magnitude in the explored density range. 

As can be seen in Fig. la and lb, the fractional abun- 
dance is a much stronger function of temperature than 
of density. This is due to the exponential dependence of 
some reaction rates on temperature. Reactions between 
neutrals that do not involve radicals or atoms often have 
considerable activation barriers. Their reactions rate is 
proportional to the Boltzmann factor ex.p{—Ea/kbT) in 
the usual Arrhenius expression, where Ea denotes the 
activation energy, ks is the Boltzmann constant and T 
the kinetic temperature of the gas. 

Selecting grid points in temperature is most important 
for the accuracy of the interpolation. In addition to 
water, sulphur bearing species also pose problems: 
Many reactions have a large activation energy and their 
fractional abundance is thus expected to be especially 
dependent on temperature. Indeed, Fig. Ic shows 
a strong dependence of H2S, CS, SO2 and atomic 
sulphur on temperature. Again, we overplot interpo- 
lated abundances onto the fully calculated results. For 
temperatures below 100 K, only a small number of 
grid points is needed, but the hot-core regime of T > 
100 K requires a much larger amount in order to keep 
deviations small. 

How can we select the grid points along the temper- 
ature axis of the grid? Two methods have been tested: 
(i) an unbiased method and (ii) hand-placed points. 
For method (i), the grid points are chosen based on the 
number of reactions with activation temperature in a 
certain temperature range. The activation temperature 
of a reaction is obtained by scaling to the activation 
energy of the very important reaction OH + H2 ^ 
H2O + H which proceeds at a gas temperature above « 
250 K (Fig. la and Stauber et al. 2006). The grid points 
are then placed to properly sample temperature ranges 
with higher activation temperatures. For method (ii), 
the grid points are placed by hand based on parameter 
studies like Fig. la. The selection of species given by 
Stauber et al. (2005) and other important molecules are 
considered at different chemical ages. 

The interpolation quality of the two methods has been 
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Fig. 1. — Fractional abundance of different species depending on physical properties (temperature, density). The arrows indicate the 
positions of the interpolation points. The red lines show the results of the interpolation in log-log space, a.) Water in a region with 
X-ray irradiation versus temperature for three different densities, b.) Density dependence of four selected species related to the ionization 
fraction, c.) Sulphur bearing species in the temperature range of the grid, d.) X-ray irradiated water at different chemical ages versus 
temperature. 



tested using the technique presented in Section 4.1. We 
find a considerably better interpolation accuracy using 
the grid with hand-placed grid points for the following 
reason: The unbiased approach does not take into 
account the relative importance of different reactions: 
e.g. the reaction between OH and H2 , despite its 
substantial influence on the chemical network since the 
key molecule water is involved. As Fig. la shows for 
this reaction, it is not suflicient to put two grid points 
near 250 K in order to obtain a good interpolation 
quality. The temporal evolution of the water abundance 
reveals another difficulty in the selection of the grid 
points (Fig. Id): For chemical ages between a few times 
10'^ and 5 x 10'' yr, the fractional abundance cannot be 
interpolated well using the same set of grid points: The 
gradient of the water abundance for 100 < T < 250 K 
and T > 250 K changes with the chemical evolution. 
Thus, the base point for the connection between the low 
and high water abundance shifts from 100 K to 250 K. 

A similar effect is also observed in Fig. la, where 
different total densities lead to different timescales of the 
water destruction. While the fractional abundance of 
water varies by approximately two orders of magnitude, 
species connected to its network may be affected even 
stronger. For example the fractional abundance of H2S 



having its main destroyer (HCO+) in common with 
H2O increases by many orders of magnitude between 
250 K and 600 K (Fig. Ic). We finally decided to 
implement 23 hand picked points at 8, 12, 17, 40, 59.9, 
60.1, 99.9, 100.1, 120, 180, 230, 250, 300, 400, 600, 800, 
1000, 1400, 1800, 2200, 2600, 3000 and 3400 K for our 
work. This range thus covers conditions from cold dark 
clouds to hot FDR like regions, heated by FUV photons. 

The selection of the grid points for the density axis is 
less critical due to the smooth abundance profile along 
this dimension. We implement one grid point per order of 
magnitude in density. The density is taken to be within 
10"' - 10^ cm""^, sufficient to model envelopes of low-mass 
and high-mass YSOs (J0rgensen et al. 2002, Maret et al. 
2002, van der Tak et al. 2000). 

2.2. FUV driven chemistry 

Reaction rates for photodissociation or ionization of 
molecules by FUV radiation can be written as (e.g. van 
Dishoeck et al. 2006) 

k= J{E)a{E)dE [s^^] . (2) 

The mean intensity of the FUV radiation field J{E) 
[cm~^ s~^ erg~^] times the cross-section for photoabsorp- 
tion and ionization (j{E) [cm^] is integrated between the 
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TABLE 1 

Parameters for the grid of chemical models, their physical range in the envelope of YSOs and the required number of 

MODELS to achieve A GOOD INTERPOLATION ACCURACY (si A FACTOR OF TWO FOR MOST SPECIES/PHYSICAL CONDITIONS). ThE 

parameters are defined AND DESCRIBED IN THE TEXT. 



Parameter 


Range 


Grid points 


straight-forward implementation: 




Density n 


lO" - 10-* cm-^ 


6 


Temperature T 


8 - 3400 K 


23 


X-ray flux Fx 


10"'^ - 10 erg cm-2 s'^ 


14 


Plasma temperature^ Tx 


10^ - 3 X 10* K 


4 


X-ray Attenuation^ Af(Htot) 


1021 _ 1024 cm-2 


6 


Cosmic-ray ion. rate Qcr 


10"" - 10-15 s-l 


5 


FUV flux Go 


- 10^ ISRF 


9 


FUV attenuation^ t 


0.1 - 10 Ay 


10 






1.7 X 10^ 


improved implementation: 






Density n 


10" - lO'' cm-3 


6 


Temperature T 


8 - 3400 K 


23 


Total ionization rate 






Ctot(Ccr, Fx,Tx,Af(Htot)) 


10"" - 10-12 s-i 


11 


a(Go,T) 


arbitrary units 


15 




arbitrary units 


5 






1.1 X 10^ 



''The models for Go=0 or Fx=0 do not depend on the attenuating 
column density and plasma temperature. The same simplification 
is possible for high values of a(Go!''"). 

average dust work function (ii^min ~ 6 eV) and the ion- 
ization energy of hydrogen (i?max = 13.6 eV). For these 
photon energies, photodissociation proceeds through line 
absorption and absorption of the continuum by dust. 
Line absorption occurs by other species or by the con- 
sidered species itself (so called self-shielding). In prin- 
ciple, the exact shape of the spectra is thus required to 
calculate the chemical rates. As an approximation, van 
Dishoeck (1988) used the spectral shape of the interstel- 
lar radiation field (ISRF, e.g. Habing 1968 or Draine 
1978) and fitted rate coefficients to the equation 

fc = Go-C-exp(-7.T) [s-'] , (3) 

where Gq is a scaling factor to the interstellar radiation 
field and t = Ai,/1.086, Ay being the extinction at 
visual wavelengths due to dust. It is calculated from the 
total hydrogen column density by the conversion factor 
7V(Htot) = 2iV(H2) + iV(H) « 1.87 x lO^i • cm-^ 
(Bohlin et al. 1978 for a dust reddening of Ry « 3.1). 
The values of the constants C and 7 in Eq. (3) depend 
on the spectrum, used to for the fitting. For hot FUV 
sources with T > 20 000 K, it is safe to use fits to 
the interstellar radiation field, dominated by O and B 
stars. Go is then scaled with respect to the average 
interstellar flux of 1.6 x 10"^ erg cm^^ s^^ at photon 
energies between 6 and 13.6 eV. For colder sources 
such as accretion hot spots, the intensity of the spectra 
decreases at lower wavelengths, and C and 7 should be 
adjusted (see van Dishoeck et al. 2008). The value of 7 
depends furthermore on the adopted dust model as well 
as on the range of Ay used for the fitting. 

Stauber et al. (2004) studied the influence of FUV 
radiation on the chemistry in envelopes of YSOs and 
found 00+ to be a good tracer in high temperature 
regions. Observations of that molecule by Stauber et al. 
(2007) revealed a surprisingly high abundance which 
cannot be explained in terms of their spherical models. 
The second part of this series of papers will thus present 
a detailed chemical model in 2D using the interpolation 



approach introduced in this work. The dependence of 
the fractional abundance of C0+ on the radiation field 
Go and the attenuation r are presented in Fig. 2. The 
fractional abundance peaks at radiation strengths of 10^ 
to IC^ times the ISRF for low extinction (r w 0.1). Self- 
shielding can lead to a high II2 abundance at extinctions 
r < 0.1. We restrict the range of our chemical grid 
to values r > 0.1, since for the planned applications, 
the uncertainty in the geometry can easily account for 
errors of r on the order of 0.1. At high extinction, with 
r > 20, the influence of the FUV radiation is negligible 
even for field strengths of 10^ times the ISRF. In 
absence of any attenuation, this FUV field corresponds 
to an early B star with a temperature of 30 000 K and a 
luminosity of 2 x 10^ Lq at a distance of about 1 000 AU. 




10° 1(1" 10' 10" 

Go |ISRF| 



Fig. 2. — Fractional abundance of 00+ for different FUV fluxes 
(Go) and FUV attenuation (r). The solid contour lines and grey 
scale give the fractional abundance of the molecule. The gas den- 
sity was chosen to be 10^ cm— ^, while the temperature is fixed at 
550 K. The curvilinear coordinate system (a, /3) given in dotted 
line is used for the interpolation in the (Go,t) - plane. The grid 
points are given by grey /red dots. 

An implementation of the chemical grid using r 
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and Go as interpolation axes would lack sufficient 
interpolation accuracy. Another approach was therefore 
developed: Contour plots of FUV sensitive molecules 
(e.g. C+, C, CO and hydrocarbons C^Hy) similar to 
Fig. 2 are calculated for different temperature regimes. 
The contour lines are then used to fit a curvilinear 
coordinate system as indicated by (red/grey) dots in 
Fig. 2. A unique function relates the physical units of r 
and Go to arbitrary units of a coordinate system denoted 
by a and (3. In this way, the interpolation quality can 
be greatly improved while the number of grid points is 
kept constant (Table 1). One grid point is placed at 
very high extinction and no FUV irradiation to rep- 
resent a model without any influence of FUV irradiation. 



2.3. X-ray driven chemistry 

Rates for the direct photoionization by X-ray photons 
can be evaluated in a similar way as for FUV radiation 
using Eq. (2). The local intensity of the X-ray emission 
of a thermal plasma is approximated by 



J{E,r) 



47rr2 ' E 



■exp 



E 



•exp ( {E) ■ Ar(Htot)) 

(4) 



where E is the photon energy, r the distance to the 
X-ray source, the temperature of the X-ray emitting 
plasma and A/'(Htot) the attenuating column density. 
The normalization factor J\f is evaluated from the 
luminosity of the X-ray source Lj. [erg s^-'^] using 
Lx = Airr"^ j I^rLa.%t.{E,r)EdE, where /unatt. denotes 
Eq. (4) without attenuation (A/'(Htot) = 0). The 
photoionization cross-section CTphoto is obtained by 
summing up the contributions of all species including 
heavy elements in the solid phase (Stauber et al. 2005). 
For a photon energy above 10 keV, inelastic compton 
scattering (described by CTcompton) governs the total 
X-ray cross-section {a^ot = Cphoto + CTcompton)- The 
influence of elastic compton scattering and line emission 
in the X-ray spectra to the chemical composition can be 
neglected (Stauber et al. 2005). 

Unlike FUV radiation, secondary processes are more 
important for the chemistry than direct photoionization 
(e.g. Maloney et al. 1996, Stauber et al. 2005): Fast pho- 
toelectrons and Auger electrons can ionize other species 
very efficiently. The rate of ionization per hydrogen 
molecule due to the impact of secondary electrons can 
be written as 



J{E,r)atotAE) 



E 



W{E)x{R2) 



dE 



(5) 



where x{H.2) « 0.5 denotes the fractional abundance of 
H2 with respect to the total hydrogen density. The mean 
energy per ion pair is given by W{E). Like cosmic rays, 
secondary electrons can excite H2, hydrogen and helium 
atoms. The electronically excited states decay back to 
the ground state by emitting FUV photons. Mostly 
photons of the Lyman- Werner bands of H2 contribute 
to this internally generated FUV field. This field can 
also photodissociate and photoionize other species and 
hence influence the chemistry. 



In a straight-forward implementation, three different 
parameters are required to parameterize the impact of 
X-rays on the chemistry: (i) The X-ray flux, as defined 
by Fx = Lx/^irr'^, with the X-ray luminosity Lx and 
the distance to the X-ray source r. {ii) The attenuating 
column density iV(IItot) and (m) the temperature Tx 
of the X-ray emitting plasma. The relevant range for 
each grid dimension is given in Table 1. It is chosen 
to cover the range of high-mass and low-mass (class 
0/1) sources as modeled by Stauber et al. (2005, 2006). 
The required number of 14 grid points in the dimension 
of the X-ray flux can be explained referring to the 
example of the very strong dependence of the water 
fractional abundance on this parameter. In Fig. 5, the 
H2O fractional abundance for a gas temperature of 120 
K is shown: At an X-ray flux of 10~^ erg s~^ cm"^, the 
abundance drops by about two orders of magnitude. As 
discussed by Stauber et al. (2006), this drop depends 
on X-ray irradiation and chemical age. Thus many 
grid points along the dimension of Fx are required to 
sufficiently sample this drop. 

Together with five points along the dimension for the 
oosmic-ray ionization rate ^cr [s~^], we would obtain 
1.7 X 10'' grid points for a straight-forward implemen- 
tation (Table 1). This large number can be reduced by a 
factor of 100 using the following approximation. Instead 
of the X-ray model described in Stauber et al. (2005), we 
make use of the fact that direct photoionization processes 
due to X-rays can be neglected for the chemical abun- 
dances: Instead of the X-ray model, the H2 ionization 
rate is calculated by Eq. (5) and the parameter for 
the cosmic-ray ionization rate is increased by this rate. 
A virtually enhanced cosmic-ray ionization rate thus acts 
as a proxy for the total ionization rate: 

Ctot = Ccr = C + Ch. {Fx,N{W), Tx) (6) 

The "standard" cosmic-ray ionization rate 
C^j. K, 5.6 X 10^^^ s^^ and the ionization rate of H2 due 
to the impact of secondary electrons (-^xj -^(H), T^,) 
form together the total ionization rate, which enters 
the chemical model as effectively increased cosmic-ray 
ionization rate. In this vein, Doty et al. (2004) used an 
increased cosmic-ray ionization rate as a proxy for an 
additional, internal source of ionization. Our approach 
treats the photodissociation reactions due to internally 
produced FUV photons exactly in the same way as the 
X-ray model by Stauber et al. (2005): Rates by Gredel 
et al. (1989) for cosmic ray induced FUV reactions are 
implemented therein with the same total ionization rate. 

This approach provides a simple recipe to include X- 
ray induced reactions with an arbitrary X-ray spectrum 
into chemical models: The average energy per ion pair 
W{E) and X-ray cross- sections (iphoto and CTcompton are 
given in Appendix A. Using Eq. (5) the total ionization 
rate C,}i^{Fx,N{]l),Tx) can be calculated. Pre-calculated 
ionization rates for a thermal X-ray spectrum (Eq. 4) 
are given in the same appendix. 



2.4. Total ionization rate vs. X-rays 

The approach to use a total ionization rate instead of 
the previous X-ray model is tested with a spherical model 



Chemical Modelling of Young Stellar Objects, I. Method and Benchmarks 



7 



of AFGL 2591 as given in Stauber et al. (2005). In this 
model, temperature and density vary with distance to 
the central source. The density structure is given in the 
form of a power-law in radius (van der Tak et al. 1999) 
and the temperature structure is taken from the detailed 
thermal balance calculations of Doty et al. (2002). The 
temperature and density distribution are given in Fig. 
3. The model is devided into 30 shells of approximately 
constant column density in radial direction. In this way, 
more shells are placed in the warm, chemically active 
part of the envelope. To resolve the water evaporation 
properly, assumed to be at 100 K, an additional shell is 
set to the corresponding position. Fewer positions are 
needed in the outer, colder and thus chemically less ac- 
tive region. Calculations with more shells show insignifi- 
cant differences in observable quantities derived from the 
models. More points may be necessary for other temper- 
ature/density profiles as presented in Fig. 3. When the 
interpolation method for chemical abundances is used to 
construct such a spherical model, the number of points 
can be easily increased due to the gain in speed using 
that approach. For the implementation, we suggest to 
bisect the points until the derived observable quantities 
converge. 




L — ■ — 

Radial Position [cm] 



Fig. 3. — Density and temperature in the spherical model of 
AFGL 2591. Vertical lines on top of the figure indicate the position 
of the shells used for the calculation. 

Figure 4 shows the radial dependence of the abun- 
dances for the main molecules CO, II2O, CO2 as well as 
for those 10 molecules predicted to be the best X-ray 
tracers by Stauber et al. (2005) (HCN, HNC, H2S, CS, 
CN, SO, HCO+, HCS+, H+ and N2H+). The ionized 
hydrides CH+, SH+ and NII+, discussed later in this 
section, are given in addition. The X-ray models are 
indicated by black lines. Results from the model with 
a total ionization rate ( given by Eq. (6) are presented 
by thick (grey/red) lines. A model without protostellar 
X-ray radiation and models with X-ray luminosity Lx 
= lO'^", 10"^^ and lO'^^ erg s~^ are shown. Protostellar 
FUV radiation is included in the same way as Stauber 
et al. (2005) assuming a flux of Go = 10 at 200 AU. The 
agreement between the model including X-rays and with 
a total ionization rate is very good for the shown set of 
species, with deviations less than 25% in the fractional 
abundance. 



This very good agreement between X-ray models and 
models with a total ionization rate has observational 
consequences: In order to use molecular lines as tracers 
for protostellar X-ray radiation, the effects of cosmic-ray 
ionization and X-ray induced ionizations have to be 
disentangled by spatial or excitation information on the 
abundance. This can be obtained from high-J lines 
with a high critical density (e.g. observed by the up- 
coming HIFI spectrometer onboard the Herschel Space 
Observatory) or by comparing visibility amplitudes from 
high angular resolution interferometer observations (e.g. 
Benz et al. 2007). 

What are the limits of this approach? The X-ray flux 
Fx = Lx/4:TTr'^ at a distance of 20 AU from a source 
with an X-ray luminosity of 10^^ erg s~^ is about 90 
erg s~^ cm~^. Meijerink & Spaans (2005) use similar 
field strengths when modeling X-ray dominated regions 
(XDRs) in galactic nuclei. Thus we explored the validity 
of the ionization rate approximation up to X-ray fluxes 
of 100 erg s~^ cm~^. All molecules in the chemical 
network have been tested for signiflcant deviations 
between the two approaches in the range from to 100 
erg s~^ cm~^. Excellent agreement as in the example 
of water (Fig. 5g) was found for most species. For our 
application in envelopes of YSOs, the ionization rate 
approach thus can be safely used. 

A small set of species revealed discrepancies between 
the two approaches under certain conditions. They 
are presented in Fig. 5. There are different reasons 
for deviation. They are discussed in the following 
paragraphs. 

The negatively charged species H~, OH^ and CN^ 
significantly deviate since the X-ray approach does not 
include a path for the formation of H" by the ionization 
of II2 (Fig. 5a). The UMIST database hsts the reaction 
II2 + c.r. H+ + for this process, resulting in a 
much higher abundance of this species. Subsequently, 
the abundances of OH^ and CN~ also increase. The 
X-ray approach only includes formation by radiative 
association (H + ^ -I- 7). Thus our total 
ionization rate approach to is an improvement. 

At low X-ray fluxes the abundance of 0+ is higher 
in the X-ray model than in the total ionization rate 
approach (Fig. 5b). In the X-ray approach, direct 
photoionization CO -f 7x ^ C+ + 0+ -I- 2e^ can 
compete with the relatively inefficient formation of O^ 
through the reaction of He+ with CO2, OCS or SO. The 
fractional abundance of 0+ thus increases due to the 
presence X-ray photons. However, it remains below a 
few times 10~^^, insignificant for our applications (Sect. 
4.1). At higher X-ray fluxes, 0+ is formed in the charge 
exchange reaction of 11+ with O, and the two approaches 
agree well. 

For the hydride ions SII+ and CH+ at low density 
and high X-ray flux, the X-ray approach predicts higher 
abundances than the ionization rate approach. This 
discrepancy can be explained by doubly ionized atoms 
added by Stauber et al. (2005) to the chemical network, 
but not included in the UMIST database. Thus doubly 
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Fig. 4. — Abundances in a splierical AFGL 2591 model calculated by the X-ray model of Stauber et al. ("X-ray", black lines) compared 
to a model with a total ionization rate given by Eq. (6) grey/red lines). Four different models are shown for X-ray luminosities of 0, 

10^", 10^1 and 10^2 erg s"!. 

ionized atoms lack in the total ionization rate approach. 
At high densities, doubly ionized atoms are not impor- 
tant because X -I- H;^ XH+ + H2 dominates XH+ 
formation. Abel et al. (2008) found the branching ratio 
of S++ + H2 to be important for the efficiency of SH+ 
production through doubly ionized atoms. Stauber et al. 
(2005) only included X++ + B2 ^ X+ + B+ with rate 
coefficients taken from Yan (1997) (fc < 10"^"^ cm^ s~^). 
However, in the later version of their code used here. 
They implemented the formation of hydride ions as the 
main product with a rate coefficient of 10"^ cm'^ s~^ 
(Yan 1997), similar to the value of Abel et al. (2008) for 
SH+, but about an order of magnitude higher for CH+. 



Doubly ionized atoms become important for low den- 
sities and high irradiation. For a given X-ray flux, the 
number of ionization processes per volume is constant, 
independent of the total gas density. The electron 
density at a high X-ray flux of 10^ erg cm~^ s~^ varies 
only by a factor of about two between a density of 10^ 
and 10* cm~^. The amount of 0+ is large, accounting 
for a short destruction time-scale of H2 in the reaction 
H2 + 0+ ^ 0H+ + H. H2 is however the precursor of 
which is subsequently also reduced (Fig. 5h). In 
this regime, the main production of SH+ and CH+ is 
through X++ + H2 ^ XH+ + H+. 
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Fig. 5. — Abundance versus X-ray luminosity of species with significant deviations between the X-ray model [black lines) compared to 
the results of the model with enhanced cosmic-ray ionization rate (grey /red lines). The temperature is indicated on the top. Models for 
different total hydrogen densities (10^, 10® and 10* cm~^) are given to point out the largest effects. Water exhibiting little deviation is 
shown as an example for comparison. 

peratm'e T, the density n, two parameters for the FUV 
flux and the total ionization rate. The implemented grid 
thus consists of five dimensions and a total of about 
f.f X fO^ grid points in the improved implementation. 
The range and number of points for each dimension are 
given in Table 1 . For the interpolation of the abundance 
at a specific physical condition A = (A^, A^, A'^, A^, A^) 
a multidimensional interpolation in logarithmic space 
is used: The neighboring interpolation points of A are 
found in each dimension d = {1,2,. ..,5}, such that 
Af < A'^ < A^. Then, the abundances at the interpo- 
lation points, x (AJ, A|, A|, Af , A^) are read out of the 
grid for all combinations of k,l,m = {a, b}. The in- 
terpolated abundance x is obtained from 

logio i^) = "fc "™ 

i,j ,k,l.'m—{a.h} 

The weights af along each dimension are defined by 
a'^ — 1 — /3 and = f3. The position within the hyper- 
cube, normalized to the interval between and 1 is given 
by /3 = [\og,,{X^) - logio(Af )] / [logio(A^) - logM] . 

The most time-consuming step of the interpolation is 
the calculation of the total ionization rate by the integral 
of Eq. (5) and the inversion of the coordinate transfor- 
mation (a,/3) {Go,t). Both functions are given in 
tabulated form in the interpolation routine. This allows 
to obtain approximately 10 000 fractional abundances 
per second on a standard personal computer. Using the 
full chemical model, the calculation of the same number 
of abundance evolutions would require about 80 hours of 
CPU time. 

3. CHEMICAL MODEL 



NH+ enhancement in the X-ray approach is due to 
another effect. The charge exchange reactions of oxygen 
and nitrogen with H2 have rate coefficients of 3 x 10^^^ 
cm^ and 10^^ cm'^ s~^, much faster than for sulphur 
and carbon. The chemical network of Stauber et al. 
(2005) follows Yan (1997) and thus assumes the forma- 
tion of 0H+ and NH+ through doubly ionized species to 
be negligible. The slightly enhanced abundance of NH+ 
in the X-ray approach compared to the total ionization 
rate approximation is explained by the larger amount of 
N"*" at low density, since its production by the reaction 
of N++ + H ^ H+ -I- N+ is faster than the reaction 
between N++ and H2. 

In conclusion, the total ionization approach should not 
be used for low densities where the X-ray fluxes exceeds 
10 erg s-i X-ray flux. For SH+ and CH+, the X-ray 
flux limit is lower (Fig. 5). We note that the X-ray 
approach is also incomplete in this regime as e.g. the 
important vibrationally excited H2 with a high energy 
deposition of X-rays per density Hx/n > 10~^^ erg cm'^ 
s""'^ (Yan 1997) is not included. Nevertheless, the total 
ionization rate approach is well applicable for models of 
YSO envelopes, since density models by J0rgensen et al. 
(2002) or van der Tak et al. (2000) predict a gas density 
larger than 10^ cm~'^ close to the protostar where X-ray 
irradiation is significant (i.e. exceeding than the cosmic 
ray effects). Indeed, the derived abundances of SII+, 
CH+ and NII+ in Fig. 4, show no difference between the 
two approaches. For our application, the total ionization 
approach is not only more elegant, but adequate. 

2.5. Multidimensional interpolation 

The relevant physical parameters for the molecu- 
lar/atomic fractional abundance evolution are the tem- 
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For our work, we start from the chemical model 
introduced by Doty et al. (2002, 2004) and Stiiuber 
et al. (2004, 2005). Thus we describe only changes 
to their model in this section: The chemical network 
has been updated from the UMIST 97 (Millar et al. 
1997a) to the UMIST 06 database for astrochemistry 
(Woodall ct al. 2007). Wc implement the standard 
UMIST 06 database without the dipolc-cnhanccd rates. 
For technical reasons, fluorine bearing species have 
not been included. Appendix B gives details of the 
implementation of the new chemical network. The 
rates for cosmic ray and X-ray induced FUV reactions 
implemented in Stauber et al. (2005) were a factor of 2 
too high due to an error in the previously used database 
(cf. comment by S. Doty in Woodall et al. 2007) and are 
now at their correct values. Charge exchange on PAH or 
small grains can be important for the ionization balance 
(Wakelam & Hcrbst 2008, Maloney et al. 1996). FUV 
irradiation enhances the number of positively charged 
grains. The rates for this process were accidentally 
divided by the total density in the implementation 
of Stiiuber ct al. (2005) and are now at the correct 
value. The influence on the chemical abundance in their 
applications is however small. Grain surface reactions 
are not taken into account except for the formation of 
H2 , where the rate given in Draine & Bertoldi (1996) is 
used. 

We have extended the chemical model with self- 
shielding of CO and H2 using the shielding factors 
given in Lee et al. (1996) and Draine & Bertoldi (1996), 
respectively. They depend on the column densities 
A''(Il2) and N{CO) between the FUV source and the 
modeled parcel of gas. Since we do not include any 
"non-local" parameter as dimension in the chemical 
grid, an approximation is needed: Prior to the cal- 
culation of a model, the optical depth-dependence 
of the II2 abundance at the given physical conditions 
(density, temperature and FUV irradiation) is calculated 
for a simple steady state-model considering only the 
II2 photodissociation and formation on dust with a fixed 
density and temperature. The column density can be 
read out of this toy model and shows good agreement 
with the examples given in Draine & Bertoldi (1996). 
As a rough approximation of the CO column density, we 
assume a fixed ratio of 7V(CO)/Ar(H2) = 2 x 10"^. 

In the papers of Doty et al. /Stauber et al., the system 
of stiff ordinary differential equations (Eq. 1) is solved 
using the DDRIV3 algorithm'^. This solver is numerically 
unstable for chemical networks involving reactions with 
short timescales, e.g. evaporation from dust or photodis- 
sociation at high values of Gq. The DVODE solver^ proved 
to be more robust and much faster in a large range of 
physical parameters. Similar characteristics in compar- 
ing the two algorithms were found by Nejad (2005). The 
equations for the conservation of elements revealed a bet- 
ter accuracy of the DVODE solver. A single run of the 
chemical model on a standard personal computer takes 
about 30 s, and the calculation of the whole grid thus 
about 900 hours of CPU time. Distribution on several 

^ http://www.netlib.org/slatec 
http://www.netlib.org/ode 



CPUs is easily possible and allows to build various chem- 
ical grids e.g. for changed chemical networks or difi^erent 
initial conditions in relatively short time. 

3.1. Initial conditions 

In order to solve the first order differential rate 
equations (Eqs. 1), initial conditions for the abundances 
have to be assumed. If the chemical evolution were 
traced starting from diffuse cloud conditions, a purely 
atomic composition would have to be assumed and the 
physical conditions such as FUV irradiation or density 
would have to be changed during the evolution (e.g. 
Lintott & Rawlings 2006). Here we follow the approach 
of Doty/Stauber et al., who started at dark cloud 
conditions with a molecular composition (Table C2 in 
the appendix). In this way, uncertainties of the physical 
evolution during this first phase and/or reactions on 
dust grains less affect the chemical composition. The 
evaporation of species frozen out on ice is not taken ex- 
plicitly into account with evaporation-type reactions but 
approximated with different sets of initial abundances 
as in the models of Doty/Stauber et al. Evaporation 
and freeze-out reactions are implemented in the model. 
However they slow down the calculation significantly 
due to short timescales at high temperature and are not 
activated for this work. 

3.2. Sulphur bearing species 

Various sulphur bearing molecules are predicted by 
Stauber et al. (2005) to be good tracers of X-ray ra- 
diation. Despite the large uncertainty in many of the 
reaction rates involving sulphur bearing species (Wake- 
lam et al. 2004), these molecules might be important 
for applications of the grid. To better constrain the ini- 
tial molecular conditions, we use the multitude of ob- 
served sulphur bearing species foimd toward AFGL 2591 
by van der Tak et al. (2003), van der Tak et al. (1999) 
and Stauber et al. (2007). Spherical models of AFGL 
2591 (cf. Sect. 2.4) are calculated for different sets of 
initial abundances, main sulphur carriers and chemical 
ages. The abundances are used as input for a full non- 
LTE radiative transfer calculation using the RATRAN 
code (Hogerheijde & van der Tak 2000) to model line 
fiuxes. A test is carried out to compare the modeled 
line fluxes with the observations. For this work we adopt 
S (T < 100 K) and H2S (T > lOOK) to be the main initial 
siilphur carriers. The abundance for the cold part follows 
Aikawa et al. (2008) to be 9.1 x 10"* relative to the to- 
tal hydrogen density. For the hot core, the abundance 
of SO is chosen to reproduce the jump in the abundance 
between the cold and hot part as observed by Benz et al. 
(2007). The adopted sulphur abundance is about the 
solar abundance (cf. Snow & Witt 1996 and Asplund 
et al. 2005) indicating no or only minor sulpur depletion 
on dust grains in the hot core. Goicoechea et al. (2006) 
modeled a PDR with a relatively low FUV irradiation 
of X = 60 ISRF (Draine) and reported a sulphur abim- 
dance of about a factor of 4 lower than the value adopted 
in this work to reproduce their observations. It may be 
caused by the massive impact of FUV irradiation in the 
vicinity of AFGL 2591, raising the temperature severely 
and resulting in an even higher sulphur abundance in the 
gas phase. 
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4. BENCHMARKS 

In this section, the accuracy of the interpolation 
method is verified. For two realistic problems - possible 
applications of the chemical grid - the interpolated abun- 
dances arc compared to abundances calculated using the 
chemical model. The results are considered satisfactory 
if they agree within a factor of two of each other for the 
reasons noted in Sect. 2. 

4.1. A spherical model of AFGL 2591 

A first test of the chemical grid is performed using 
the spherical model of AFGL 2591 as introduced in 
Sect. 2.4. Figure 6 shows the radial dependence of the 
abundances for the main molecules CO, H2O, CO2 as 
well as for those 10 molecules predicted to be the best 
X-ray tracers by Stauber et al. (2005) (HON, HNC, 
H2S, CS, ON, SO, HCO+, HCS+, H+ and N2H+). 
C0+ which will bo addressed in the next paper of this 
scries and the two molecules, CbH and HCNH+, having 
the largest deviation arc shown in addition. Solid lines 
give the results of the chemical model, while interpolated 
abundances are shown by dashed lines. Models with no 
protostcUar X-ray radiation and an X-ray luminosity of 
Lx — lO''^ erg s^^ (Stauber et al. 2005) are given in thin 
and thick lines, respectively. We assume a chemical age 
of 5 X 10^ yr. The shaded region indicates a range of a 
factor of two compared to the fully calculated model and 
marks the goal of our interpolation approach. Indeed, 
most molecules comply with the aimed accuracy for the 
models with and without X-rays. 

To provide an unbiased check for a larger set of 
molecules, a statistical approach is used. The goal of 
chemical modeling is to obtain characteristics which can 
be compared to observations. In the following, we check 
the column density of all molecules in the chemical model 
for agreement between grid and full calculation. The ra- 
dial column density, iV^adiai = / n{r)dr, with the radial 
distance r, gives a measure for molecules observed in ab- 
sorption. For molecules observed in emission, the beam 
averaged column density is more appropriate. It is de- 
fined by 



A^beam — 



J J n{z,p)G{p) 2Trpdpdz 
/ G{p) 2'Kpdp 



(7) 



with z being the coordinate along the line of sight and p 
the impact parameter perpendicular to z. G{p) is a beam 
response function (e.g. a gaussian). This approach using 
column densities does not incorporate the molecular 
excitation. Thus it cannot be used to compare models 
directly with observations. Computationally much more 
demanding radiative transfer calculations would have 
to be carried out instead. For our benchmark purposes 
however, column densities are sufiicient. 

Four different spherical models of AFGL 2591 with- 
out protostcUar X-rays and with X-ray luminosity of 
Lx = 10''°,10''^ and lO"*^ erg s^^ are compared. Three 
difi'erent chemical ages of 3 x 10^, 5 x lO'' and 3 x 10^ 
yr are considered to trace possible deviations due to the 
temporal evolution of the fractional abundances. For 
these 12 models, the radial (Aj-adiai) and beam averaged 
(A^beam) column density are obtained for all molecules 



in the chemical network. We assume a 14" beam 
corresponding approximately to a JCMT or Herschel 
beam and adopt a distance of 1 kpc to the source. To 
compare grid results to fully calculated models, the 
factorial deviation Y = max (AfgHd/Afuiii A^fuu/A^grid) is 
introduced. The fraction of molecules having a factorial 
deviation larger than a certain value is shown in Fig. 
10. Only molecules with a column density larger than 
10^-^ cm~^ are considered for the following reason: 
Assuming optically thin radiation, the line strengths 
corresponding to this column density for a molecule 
at maximal excitation, a typical Einstein-A coefficient 
of 10""^ s~^ and a frequency of 350 GHz is about 40 
mK km s~^, a lower limit for observations with current 
and upcoming facilities. Compared to the H2 column 
density of order 10^^ cm~^, it corresponds to a fractional 
abundance of about 10^^^ which marks approximately 
the possible limit of a detection within an observation 
time of a few hours. 

Table Dl in Appendix D lists the molecules with the 

largest deviations. The largest disagreement is found 
in the radial column density of OH^ with a factorial 
deviation up to a factor of 9.6. This deviation is 
explained by the different paths for H~ formation in the 
two models (Sect. 2.4). It is thus an improvement of the 
total ionization approach. A problematic disagreement 
larger than a factor of 5 is found however for IICNII+ 
(r = 5.4) and CeH (Y = 5.2). Figure 6 shows the radial 
dependence of the fractional abundances for these two 
species at the chemical age and X-ray luminosity where 
the largest deviation occurs. 

How can such large deviations be explained? In the 
case of CeH at low X-ray fluxes, it is surprisingly an 
effect of the ionization rate. Figure 7 presents the 
fractional abundance of CgH vs. total ionization rate 
obtained from the chemical model and interpolated from 
the chemical grid. Three different chemical ages are 
given. For 5 x 10'' yr, we find a deviation of about a 
factor of 5 due to the interpolation over a sharp peak 
in the fractional abundance (bold vertical bar), which 
is not resolved by the grid. This peak shifts to higher 
ionization rates with time. A proper sampling would 
thus require a large number of grid points to cover all 
chemical ages. At ionization rates below this peak, 
CgH is formed by the reaction C + C5H2 CgH. C 
stems from the X-ray or cosmic ray induced dissociation 
of CO. At an older chemical age (or higher ionization 
rate), the "late-phase" molecule SO2 is photodissociated 
to SO and O. Reaction of CgH -|- O then decreases the 
fractional abundance of CgH. At very high ionization 
rate, the bulk of CeH is produced by the electron 
recombination of C7H+. 

Most contribution to the radial column density of 
HCNH+ is from the innermost part of the model. In 
this region, HCNH"*" is mainly formed by the reaction 
of HCN with H3O+. Both molecules arc enhanced by 
X-rays and FUV irradiation at temperatures above 250 
K, where water is not destroyed due to X-ray irradiation 
(Stauber et al. 2006). An X-ray luminosity higher than 
10^° erg is sufficient to dominate the effect of FUV 
enhancement. At lower X-ray luminosities however, the 
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Fig. 6. — Abundances of major species, X-ray/UV tracers and problem cases in a spherical model of AFGL 2591. The plots show a 
compajrison between results calculated by the chemical model (solid line) and interpolated from the grid of chemical models (dashed line). 
The grey region shows a deviation of a factor of 2 from the chemical model. Note the different vertical scales of individual panels. The thick 
lines corresponds to a model with a stellar X-ray luminosity of 10^^ erg s~^ and the thin line to a model without X-rays. All abundances 

are given for a chemical age of 5 X 10^ yr, except for HCNH+ (3 X 10^ yr). 



FUV flux dominates the HCNH+ abundance. Due to 
the particular shape of the contour lines of its fractional 

abundance depending on Gq and r, this molecule is 
difficult to interpolate. Figure 8 shows the fractional 
abundance of HCNH+ in the (Go, T)-plane along with 
the coordinate system used for the gridding and interpo- 
lation (Sect. 2.2). Physical conditions of the innermost 
region of the AFGL 2591 model and a chemical age of 
3 X 10^ yr were chosen for this figure. 

In general, the agreement between the grid results 
and the fully calculated abundances is excellent. The 

mean deviation of (Yboam) = 1-12 for the beam aver- 
aged column densities and (^radial) = 1-28 for the radial 
column densities indicate agreement for the majority of 
the species. A total number of 1661 (A^radiai) and 1563 



(A^beam) columu densities are considered for the means, 
respectively. The slightly higher mean deviation of the 
radial column density is explained by the larger weight 
of the chemically active hot-core region. The mean de- 
viation does not vary by more than 0.04 with chemical 
age. 

4.2. FUV driven chemistry 

In the spherically symmetric models of AFGL 2591, 
protostellar FUV radiation cannot penetrate further into 
the envelope than a few hundred AUs (Stauber et al. 
2004). This is however not the case when an outflow 
cavity allows FUV photons to escape and irradiate a 
larger volume of gas. This 2D situation will be addressed 
in the second and third parts of this series of papers. 
Here we carry out the following benchmark test. A 
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Fig. 7. — Illustration of the CeH problem: Fractional abundance 
of CeH vs. ionization rate calculated by the full chemical model 
(black lines) and from the grid of chemical models (grey/red lines). 
Three different chemical ages are given. The grid points are indi- 
cated at the top of the figure. The cosmic-ray ionization rate of 
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s is given by the dotted line. 
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Fig. 8.— Illustration of the HCNH+ probl em: Fractional abun- 
dance of HCNH"'' (grey scale and contours) at a chemical age of 
3 X 10^ yr and physical conditions corresponding to the innermost 
part of the AFGL 2591 model with Lx = 10^" erg s"!. The dotted 
grey /red lines indicate the coordinate system used for the interpo- 
lation and the solid points correspond to grid points. 

(plane parallel) region with a fixed temperature of 80, 
100 and 550 K and a density of 10^ cm~'^ is considered. 
The incident radiation field is assumed to be 1, 10 10^, 
10^ and 10^ times the ISRF and no geometrical dilution 
(i.e. Go oc r^^) is taken into account. The dust column 
density however attenuates the FUV radiation. It is 
given by the optical depth r using the conversion factor 
introduced in Sect. 2.2. Figure 9 shows the fractional 
abundance vs. the optical depth r. The main molecules 
CO, H2O, CO2, several PDR related species and the 
molecules showing the largest deviations are selected 
for Fig. 9. For most species only the results for a 
temperature of 550 K are shown, since FUV irradiation 
can considerably heat the gas through the photoelectric 
effect on dust grains. For reasons of clarity not all 
fractional abundances are given in this plot. 

A value reflecting observable quantities is needed for 



the statistics comparing interpolated to fully calculated 
abundances (Fig. 10). We use the column density of a 
molecule i between r = 0.1 and r = 15 calculated by 



N, 1.87 X 10^ 



cm 



(8) 



0.1 



where Xi(r) is the fractional abundance relative to the 
total gas density of the species i depending on the 
optical depth r. All molecules in the chemical network 
are considered at five different values of Go and three 
temperatures. Only column densities larger than 10^^ 
cm~^ are selected following the argument given in Sect. 
4.1. A total of 2096 column densities are taken into 
account. The mean deviation for this benchmark is 
found to be (Ifuv) = 1-35 and thus higher than the 
values in the previous section, but well within the goal 
of a factor of two. 

Table Dl in Appendix D lists species, for which large 
factorial deviations in the column density between the 
grid and the fully calculated model were found. HSO"*" 
shows the largest deviation amounting to a factor of 
Y = 5.6 for a high FUV flux of Go = 10*^ times the 
ISRF at a temperature of 100 K. The deviations can be 
explained by the sharp peak in fractional abundance at 
r « 11. It is a result of the formation through SO + 
^ HSO+ + H2 or SO + HCO+ ^ HSO+ + CO. The 
peak in the abundance is also found in the fractional 
abundance of SO due to the competition of two FUV 
related reactions: At r < 11, SO is formed by the 
photodissociation of SO2 and destroyed by the reaction 
with OH leading to SO2. At larger r however, SO 
dissociation becomes more important and its abundance 
thus decreases. This peak is not resolved by the current 
implementation of the grid as can easily be seen in a 
plot of the fractional abundance depending on Go and 
r (Fig. Dl in Appendix D). While large deviations of 
IISO+ occur at high r, the Oj fractional abundance 
disagrees at low optical depth (Fig. D2). 



What is the reason for these peaks in fractional 
abundance? Abundances in FUV irradiated re- 
gions are controlled by rate coefficients of the form 
k oc exp (—7 • r), where the product of the optical depth 
times a coefficient enters (see Eq. 3). Photodissociation 
rates thus drop exponentially with t. They become 
unimportant for the chemical abundance in a narrow 
interval of r. It should be noted however, that 7 is also 
affected by uncertainties at a level of at least 10% due to 
the dust model used for fitting the rate and the range of 
r, over which the rates have been fitted (van Dishoeck 
et al. 2008). 

To reveal areas in the (Go,t) plane where devia- 
tions occur, the mean factorial deviation for all three 
temperatures and all species is shown in Fig. 11. It 
shows that large areas agree well with (Y) < 2, while 
other regions show a larger disagreement, however still 
mostly within a factor of 4. A narrow region in the 
(Go,t) plane is found to have mean deviations larger 
than a factor of 10. It coincidences with the values of 
Go and r where a large deviation of IISO+ has been 
found in a previous paragraph. As the statistics in Fig. 
10 and the low mean deviation of the column density 
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Fig. 9. — Fractional abundances of FUV sensitive species between r = and r = 15. For this presentation, the temperature is fixed at 
the given value and the density is assumed to be 10 cm' Results for an incident FUV field of 1, 10^, 10* and lO'* ISRF are shown. The 
thick line corresponds to the result calculated by the chemical model while the thin line gives the results from the grid interpolation. The 
grey region indicates the range of a factor of two relative to the calculated model specified for agreement of the interpolation. 



show, this local deviation is however sufficiently averaged 
out when observable quantities are derived from the grid. 

5. UTILIZING THE GRID OF CHEMICAL MODELS 

The grid of chemical models is accessi- 
ble through AMOUR (Abundance Modelling 
of young stellar Objects Under protostel- 
lar Radiation) available for public use at 
http: / / www.astro.phys.ethz.ch/ chemgrid.html! . The 
page provides an online form to retrieve interpolated 
fractional abundances from the grid of chemical models. 
This form also calculates the X-ray ionization rate 
for a particular X-ray flux Fx, X-ray attenuation 
■^(Htot) and plasma temperature Tx- For large modeling 
tasks or to include the interpolation method in other 
codes, the page also offers a FORTRAN code and the 
necessary molecular data tables in binary and ASCII 
format for download. Documentation is available at the 



website which describes the input/output parameters, 

the technical implementation, the format of the data 
tables, and shows example input/output data. 

There exist a number of restrictions for the currently 
provided implementation which must be kept in mind 
when applying the method: The chemical composition in 
the molecular data table has been calculated for a specific 
chemical network (the UMIST 2006 database,Woodall 
et al. 2007) using one set of initial conditions (Table 
C2). The current implementation does not include grain 
surface reactions except for H2 formation on dust. The 
medium is assumed to be static and physical conditions 
(e.g., temperature or density) have been assumed to be 
constant over time. These are however only restrictions 
of the currently distributed molecular data tables. Fur- 
ther tables for other chemical networks, initial condi- 
tions, time-dependent parameters, etc. can be obtained 
from the authors on a collaborative basis. 
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Fig. 10. — statistics on the accuracy of the interpolation method 
for chemical abundances: The fraction of species within a factorial 
deviation given in the X-axis is shown. Both benchmarks to the 
spherical AFGL 2591 model (Sect. 4.1) and the FUV model (Sect. 
4.2) are given. For the AFGL 2591 model, the comparisons of 
radial and beam-averaged column densities are shown. 
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Fig. 11. — Mean factorial deviation Y between the interpolated 
and fully calculated abundances of all species in the chemical net- 
work (grey scale). Go is the FUV field strength in units of the 
ISRF and t the optical depth to the FUV source. The curvilin- 
ear coordinate system (a, f3) given in dotted line is used for the 
interpolation in the (Go,t) - plane. The grid points are given by 
grey /red dots. 

6. CONCLUSIONS AND OUTLOOK 

Starting from the chemical model of Doty ct al. (2002, 
2004) and Stauber et al. (2004, 2005), we have intro- 
duced a new method to simulate the chemical evolu- 
tion of YSO envelopes based on interpolation in a pre- 
calculated grid. The chemical model has been revised 
and the relevant physical parameters for the chemical 
composition of the gas are discussed. Benchmark tests 
have been carried out to verify the accuracy of the inter- 
polation method. We conclude the following: 

1. Accurate chemical modeling of a multidimensional 
envelope of YSOs is possible using a fast, grid- 
based interpolation method (Sect. 4). Comparison 
of observable quantities (beam averaged line fluxes) 
yields a mean factorial deviation between fully cal- 
culated and interpolated values of 1.35. This is 



more accurate than the uncertainties introduced 
by observations and chemical rate coefficients. Our 
method is more than five orders of magnitude faster 
than the full calculation (Sect. 2.5). 

2. The X-ray models by Stauber et al. (2005) can be 
reproduced using an enhanced cosmic-ray ioniza- 
tion rate as a proxy for X-ray irradiation (Sect. 
2.4). In the relevant physical range for the applica- 
tion in YSO envelopes, the agreement in the frac- 
tional abundance is within 25%. The implemen- 
tation of the this approach is described in detail 
in Appendix A and allows to include the effect of 
X-ray irradiation to chemical models in a simple 
way. 

3. Ionization by X-rays and cosmic rays cannot be 
easily distinguished by molecular tracers. Spatial 
information on the abundance is thus needed to 
disentangle protostellar X-ray and cosmic-ray ion- 
ization (Sect. 2.4) 

4. For the formation of CII+ in low density gas (lO'' 
cm~'^) with high X-ray irradiation (> 1 erg ~^ 
cm~^) recombination of doubly ionized carbon with 
H2 (C++ + H2 ^ CH+ + H+) is important. 

5. Increasing the initial abundance in the hot-core re- 
gion (T > 100 K) of the main sulphur carrier im- 
proves the agreement between models and obser- 
vations in a high-mass star forming region (Sect. 
3.2). 

6. Exploring the parameter space, we find regions 
with high gradients in molecular abundances. This 
is where small changes in the physical parameters 
yield large variations in abundance. As an exam- 
ple, interpolation of fractional abundances of an 
FUV irradiated region is difficult because photodis- 
sociation rates depend exponentially on the optical 
depth (Sect. 4.2). We have compensated for this 
by adopting a curvilinear coordinate system and a 
high number of grid points (Sect. 2.2). 

The current grid method is limited by the assumption 
of the initial composition and fixed chemical rates. 
For different assumptions, the database needs to be 
recalculated. Using this fast chemical interpolation 
method, the radiative transfer calculation becomes 
the bottleneck in computing time to interpret data. 
Recently introduced escape probability methods, such 
as the exact method by Elitzur & Asensio Ramos (2006) 
or an approximative multidimensional code by Poelman 
& Spaans (2005) can speed up this step of the modeling. 

This paper shows the possibility of interpolation for 
chemical modeling. In future publications, we will 
demonstrate major applications in multidimensional 
chemical modeling of YSOs and fast data-fitting to inter- 
pret observations. The gain in speed allows to carry out 
parameter studies on the influence of the geometry on 
the interpretation of observations. Furthermore, it will 
be possible to apply detailed chemical models to a large 
set of sources and draw conclusions on physical proper- 
ties based on statistics. 
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APPENDIX 

A. CALCULATION OF THE IONIZATION RATE 
For the calculation of the X-ray ionization rate using the equation 

Ch. = f J JiE, rME)^,^£-^^dE , (Al) 

we provide the compton and photoionization cross-sections as well as the mean energy W{E) of an ion-pair in a 
tabulated form. The cross-sections are calculated using the X-ray model of Stauber et al. (2005) to which we refer for 
the source of the atomic and molecular constants and the exact implementation. We use the elemental composition 
given in Appendix C. Since the photoionization cross-section of a molecule is approximated by adding the cross-sections 
of the contained atoms, the total photoionization cross-section does not depend on the abundance of the molecular 
species as long as the elemental composition does not change. 



T I 1 1 I 1 1 




Energy [eV] 

Fig. Al. — Photoionization and compton cross-section versus plioton energy. The ticks at the bottom of the figure indicate the supporting 
points for the interpolation of cr given in Table Al. 



The cross-sections are shown in Fig. Al for a photon energy between 100 and 10^ eV. The photoionization cross- 
section, given by the solid line is approximated using a power-law, defined piecewise for an energy range between E'min 
and i^max- Table Al gives the cross-sections at the boundaries of each energy interval. The cross-section at an energy 

E G [i^min , -Eniax] Can then be calculated using 

CTphoto(^) » 10l°Sio('^-'»)-(l-«)+l°Sio(<^-ax)-a ^ (A2) 

with a = (\og{E) — log(i?,„iji))/(log(_Ei„ax) — log(-Emin)). The deviation of the fit to the calculated cross-section is less 
than 5% in the given energy range. Interpolation intervals are indicated by ticks at the bottom of the figure. 

The cross-section for inelastic compton scattering is dominated by H2 and H. This process does not contribute to the 
total cross-section at low energy. We therefore use the fit from Stauber et al. (2005) to the XCOM database (NIST) 
for the energy above 1 keV. With x = logiQ{E [eV]), the cross-section reads 

cTcompton = 2.869674 • 10~23 _ 2.6364914 • lO^^s . ^ (A3) 
-h7.931175 • lO-^"* • - 7.74014 ■ 10"^^ • 
[era] (i;<10keV) 
CTcompton = -2.374894 • 10"^^ -I- 1.423853 • 10"^* • x 

-1.70095 • 10"^^ • a;^ [cm^] {E > 10 keV) 
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Finally, the mean energy per ion pair W{E) is constant above 1 kcV, W{E > l[keV]) = 20.95 
approximated by 

W{E) = 23.65 - (logio(i;[eV]) - 2) • 2.7 eV 
between 100 and 10^ eV, following Dalgarno et al. (1999). 



TABLE Al 

Fitting parameters for the X-ray cross- section. a(b) means ax 10*'. 



Emin [eV] 


Emax [eV] 


f^min 


[cm'^] 


(^max 


[cm-^] 


100 


291 


6.02( 


-20) 


2.71( 


-21) 


291 


404.7 


3.03( 


-21) 


1.15( 


-21) 


404.7 


538 


1.22( 


-21) 


5.29( 


-22) 


538 


724 


7.97( 


-22) 


3.59( 


-22) 


724 


857 


3.99( 


-22) 


2.54( 


-22) 


857 


870 


2.59( 


-22) 


2.48( 


-22) 


870 


1311 


2.89( 


-22) 


9.75( 


-23) 


1311 


1846 


1.06( 


-22) 


4.13( 


-23) 


1846 


2477 


4.61( 


-23) 


2.03( 


-23) 


2477 


3203 


2.23i 


-23) 


1.09( 


-23) 


3203 


4043 


1.11( 


-23) 


5.76( 


-24) 


4043 


5996 


5.89( 


-24) 


1.89( 


-24) 


5996 


7124 


1.91( 


-24) 


1.15( 


-24) 


7124 


8348 


2.24( 


-24) 


1.45( 


-24) 


8348 


28900 


1.50( 


-24) 


4.24( 


-26) 


28900 


100000 


4.24( 


-26) 


1.04( 


-27) 



THE IONIZATION RATE OF A THERMAL SPECTRUM 

As an important example, wc calculate the ionization rates ( [s^^] for a thermal X-ray spectrum. The photo- and 
compton ionization cross-sections arc implemented in the way described before and Eq. (Al) is evaluated. Fig. A2 
shows ( depending on the attenuating column density. The "standard" cosmic-ray ionization rate of 5.6 x 10~^^ s~^ 
is given by a dotted line. The top panel gives the ionization rate for 3 different plasma temperatures at a fixed flux 
of 10~^ erg s^^ cm~^. A hotter spectrum results in a larger number of photons at high energy which can penetrate 
further into the envelope. Most photons of a cold spectrum (« 10^ K) however are quickly absorbed. For the 
bottom panel, the plasma temperature has been fixed to 7 x 10^ K. Since the X-ray intensity depends linearly on 
the flux _Fx, the ionization rates scale in the same way. The difference in the ionization rate, when the integral in 
Eq. (Al) is evaluated between 1-100 keV (solid line) and 0.1-100 keV (dashed line) is small at column densities 
larger than 10^^ cm~^, except for a small difference due to the normalization J\f (Eq. 4). Up to column densities of 
a few times 10^^ cm~^, photons with an energy below 1 keV can however contribute significantly to the ionization rate. 

In Table A2, the ionization rate for an X-ray flux of 1 erg s~^ cm~^ is given. The columns correspond to different 
plasma temperatures and the rows to different attenuating column densities. These ionization rates can be scaled 
linearly to an arbitrary value of the X-ray flux. For a point-hke X-ray source, the flux at a distance r from the source 
is obtained from the X-ray luminosity by Fx = Lx/A-Kr"^. 



TABLE A2 

The ionization rate in at Fx = 1 erg s^-'^ cm^^. a(b) means a x 10'\ The columns give values for different plasma 

temperatures [K] and the rows corresponds to C;0LUMN densities between 10-^° AND 10^^ cm~2. 



Af(Htot) 




3(6) K 


1(7) K 


3(7) K 


1(8) K 


3(8) K 


1(20) 


cm" 


-•1 


4.8(-ll) 


2.4(-ll) 


l.O(-ll) 


3.4(-12) 


1.2(-12) 


1(21) 


cm" 


-2 


4.2(-12) 


4.5(-12) 


2.5(-12) 


9.8(-13) 


3.7(-13) 


1(22) 


cm" 


-2 


6.9(-14) 


3.7(-13) 


4.0(-13) 


2.2(-13) 


1.0(-13) 


1(23) 


cm" 


-2 


6.9(-17) 


9.6(-15) 


3.9(-14) 


4.5(-14) 


3.3(-14) 


1(24) 


cm" 


-2 


6.8(-22) 


4.5(-17) 


2.0(-15) 


9.4(-15) 


1.4(-14) 


1(25) 


cm" 


-2 


1.8(-30) 


8.6(-22) 


1.7(-17) 


1.4(-15) 


6.4(-15) 
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Fig. A2. — Ionization rates assuming a thermal X-ray spectrum. The X-axis gives the attenuating column density. Top: Ionization rate for 
different plasma temperatures for an X-ray flux of 10"'^ erg s~^ cm~^. Bottom: Ionization rate for different X-ray fluxes (10~^, 10~^,0.1 
and 10 erg cm~^). Solid line: Eq. (Al) is integrated from 1 - 100 keV. Dashed line: Eq. (Al) is integrated from 0.1 - 100 keV. 
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TABLE Bl 

Reactions not extrapolated to lower temperature in the 

IMPLEMENTED NETWORK (CF. TEXT). THE RATES AT 10 AND 100 K 
ARE OBTAINED USING THE UMIST 06 DATABASE. THE NUMBER OF THE 
REACTION IN THE UMIST DATABASE IS GIVEN IN THE FIRST COLUMN. 



Nr. Reaction fc(lOK) fc(lOOK) 

[cm'"' s^^] [cuV' s^^] 



158* 


CH + CH3OH -» CH3 + H2CO 


1.8 


-7) 


2.1 


-9) 


239"^ 


CH2 + CH2 ^ C2H3 + H 




3.3 


-11) 


3.3 


-11) 


242^ 


CH2 + -» CO + H2 




8.0 


-11) 


8.0 


-11) 


243'' 


CH2 + -* HCO + H 




5.0 


-11) 


5.0 


-11) 


291'' 


CH3 + OH ^ H2CO + H2 




1.7 


-12) 


1.7 


-12) 


301^* 


CH3 + 02^ HCO + H2O 




1.7 


-12) 


1.7 


-12) 


338'' 


-1- HoCO — > CO + OH -1- H 




1.0 


-10) 


1.0 


-10) 


420'' 


NHo 4- NO ^ No -1- OH -1- H 




1.5 


-12) 


1.5 


-12) 


430"^ 


OH + CN OCN + H 




7.0 


-11) 


7.0 


-11) 


445'' 


OH -1- HOCH — > POo -1- Ho 4- 




1.8 


-11) 


1.8 


-11) 


457'' 


NH3 + CN ^ HCN + NH2 




2.0 


-9) 


1.0 


-10) 


491'' 


nivT -1- r^,, ^ or^i\T _i_ 

K^i>i -\- yJ2 — * -\- KJ 




1.6 


-9) 


3.2 


-11) 


515'' 


HCO + HCO ^ CO + CO + 


H2 


3.6 


-11) 


3.6 


-11) 


4084'' 


H + OH -» H2O + 7 




3.3 


-14) 


6.6 


-16) 


4552<^ 


CO + M-*0 + C + M 




< 1 


.0{-30) 


< 1 


.0(-30) 


4554' 


H + CH3 CH4 + M 




2.8 


-26) 


4.5 


-28) 


4555' 


H + ^ OH + M 




1.3 


-30) 


1.3 


-31) 


4558' 


H + OH ^ H2O + M 




1.3 


-28) 


1.4 


-30) 


4563' 


H + O2 -> O2H + M 




2.2 


-28) 


1.2 


-31) 


4568'' 


H2 + N ^ NH2 + M 




1.0 


-26) 


1.0 


-26) 


4571*^ 


C + C ^ C2 + M 




1.9 


197) 


1.7 


-9) 


4573=* 


C + ^ CO + M 




4.9 


67) 


9.6 


-19) 


4574=* 


C+ + CO+ + M 




4.9 


69) 


9.6 


-17) 


4575=^ 


C + 0+ ^ CO+ + M 




4.9 


69) 


9.6 


-17) 


4579=^ 


+ O2 + M 




3.2 


-9) 


9.4 


-32) 


4581' 


+ SO SO2 + M 




4.8 


-28) 


6.9 


-30) 


4582^ 


OH + OH ^ H2O2 + M 




2.2 


-2) 


1.9 


-28) 


4583=- 


O2H + O2H -> H2O2 + O2 + 


M 


2.8 


10) 


3.6 


-29) 



''Kept constant below 500 K, k{T < 500K) = fc(500K). 
l-Kept constant below 1000 K, k{T < lOOOK) = fe(lOOOK). 
■^Switched ofT below 500 K, k(T < 500K) = 0. 

''Switched off outside recommended temperature range, k{T 

['^min, '^max]) = 0. 

''Also listed in the OSU database, thus kept at recommended rate. 
'Kept constant below T^in, k{T < T^in) = K'^min)- 

B. IMPLEMENTATION OF THE UMIST 06 NETWORK 

The UMIST 06 database for astrochemistry (Woodall at al. 2007) lists an interval [Tmin, Tmax] for the recommended 
temperature range for each reaction rate. For some reactions, different rates are given for distinct temperature ranges. 
The paper lists several rates, for which no extrapolation to low temperatures of w 10 K should be done. We switch 
off those reactions for temperatures below Tmin- 

All other reaction rates are extrapolated to lower temperatures with a few exceptions: Some reactions of the collider 
type ( "CL" in Woodall et al. 2007) have a negative activation energy. To stabilize the chemical network, we keep their 
rates constant at temperatures below 500 K or 1000 K, or switch the reaction entirely off below 500 K (cf. Table Bl). 

Other reactions involve significant extrapolation, i.e., T^i^ > 300 K, and the rates do not decrease from 100 K to 10 
K. Certain reactions are thus switched off outside their temperature range as given in the database. The reaction OH 
+ CN — » OCN + H is not switched off, since the OSU database^ (Smith et al. 2004) lists this reaction with the same 
rate description. The reaction O + H2CO — > CO + OH + H was introduced in UMIST 99 and has severe consequences 
for the abundance of formaldehyde: it is mainly destroyed by this reaction and the abundance drops by about two 
orders of magnitude at steady state conditions. The line flux, modeled with RATRAN (Hogerheijde & van der Tak 
2000) of a spherical model of AFGL 2591 is short by about two orders of magnitude compared to observations by van 
der Tak et al. (1999). Thus we switch this reaction off below the recommended temperature range of 1750 K - 2575 K. 
Furthermore reactions with a rate more than an order of magnitude faster at 10 K than at 100 K are kept at constant 
below their minimum temperature, k{T < Tmin) = ^(Tmin)- 



www.physics.ohio-state.edu/~eric/research.html 
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C. INITIAL CONDITIONS 

The total elemental abundances in Table CI are used for the calculation of the photoionization cross-section, where 

also heavy elements locked into dust-grains arc taken into account. The values are taken from Yan (1997), except for 
helium, where the values assumed by Stauber et al. (2005) have been adopted. 

The initial conditions of species in the gas phase relative to the total hydrogen density (ntot = 2n(H2) + n(H)) are 
given in Table C2. If no other reference is given, we follow Doty et al. (2002, 2004) and Stauber et al. (2004, 2005) 
for the abundances. All species in the network without specification are initially set to an absolute abundance of 
10~* cm~"^ (effectively zero). 



TABLE CI 

Total elemental composition. a(b) means a x 10*'. 



Element 


ntotW/ntot 


Element 


«tot(X)/ntot 


H 


1.0 


Fe 


3.2(-5) 


He 


8.5{-2) 


Ne 


1.4(-4) 


C 


3.5(-4) 


Na 


2.1(-6) 


N 


1.0(-4) 


Mg 


4.0(-5) 


O 


5.4(-4) 


Al 


3.1(-6) 


s 


2.0(-5) 


Ar 


3.8(-6) 


p 


1.0(-8) 


Ca 


2.2(-6) 


Si 


3.5(-5) 


Cr 


4.9(-7) 


CI 


8.3(-8) 


Ni 


1.8(-6) 
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TABLE C2 
Initial conditions in the gas-phase. a(b) 

IIEAXS a X 10''. 



Species 


nGaB{X)/ntot 




Remark 




0.5 






CO 








C02 




T < 100 K 


a 






T > 100 K 




H20 




T < 100 K 


a 




7.5(-5) 


T > 100 K 







4(-5) 


T < 100 K 








T > 100 K 


a 


H2S 


_ 


T < 100 K 


a 




2(-5) 


T > 100 K 


b 


s 


9 l(-8,) 


T < 100 K 


b,c 






T > 100 K 


a 


N2 


3.5(-5) 






CH4 


5(-8) 






C2H4 


4(-8) 






C2H6 


5{-9) 






H2C0 




T < 60 K 


a 




6{-8) 


T > 60 K 


d 


CH30H 




T < 60 K 






5(-7) 


T > 60 K 


d 


He 


8.5(-2) 






He+ 


2(-10) 




e 


H 


5(-8) 




e 


H+ 


3(-10) 




e 




2(-9) 




e 


HCO+ 


3(-9) 




e 


H3O+ 


5(-10) 




e 


Grain" 


1.9(-8) 




e,f 


e~ 


7.5(-9) 




e 



''assumed to be frozen-out onto dust grains at 
cold temperatures or to be not abundant in hot 
regions. 

''see Sect. 3.2 

•^Aikawa et al. (2008) 

■^Evaporation temperature from Doty et al. 
(2004) 

"'Stauber et al. (2004, 2005) 

^Negatively charged grains (Maloney et al. 1996) 
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TABLE Dl 

List of molecules with the largest factorial 
deviation of the grid interpolation from the fully 

CALCULATED MODEL FOUND IN BENCHMARK TESTS IN SeCT. 

4. The column densities of both approaches are 
given along with the factorial deviation y. 



Species 


[cm^^] 


[cm "^l 


Y 






Sect. 4.1: AFGL 2591 model (A^radial) 


Lx 


Age 










[erg s"-'] 


[yrs] 


OH- 


1.9(10) 


1.8(11) 


9.6 


10^2 


3(3) 


OH- 


1.1(11) 


8.9(11) 


7.9 


10^1 


3(3) 


CeH 


3.2(12) 


6.2(11) 


5.2 





5(4) 


HCNH+ 


6.4(13) 


1.2(13) 


5.1 


10=^° 


3(5) 


C7H 


2.1(12) 


4.2(11) 


5.0 





5(4) 


HCOOCH3 


5.7(11) 


1.1(11) 


4.9 





5(4) 


H- 


1.3(11) 


5.9(11) 


4.5 




5(4) 


HS+ 


8.5(11) 


2.0(11) 


4.3 




o^^o J 


H- 


1.4(11) 


6.0(11) 


4.3 




3(5) 


H- 


1.5(11) 


6.2(11) 


4.2 




3(3) 


CH2CHCN 


2.2(11) 


5.2(10) 


4.2 


10^1 


5^4) 




4.9(12) 


1.2(12) 


4.0 


n 


0\0 J 


HS+ 


8.5(11) 


2.1(11) 


4.1 




5(4) 


Sect. 4.1: AFGL 2591 model (Afbeam) 


r 


Age 










[erg s ] 


[yi's] 


CH3OH 


8.3(12) 


2.7(12) 


3.0 





5(4) 


CH3OCH3 


5.7(11) 


2.0(11) 


2.8 





5(4) 


H2S 


2.2(15) 


1.0(15) 


2.2 





5(4) 


HC7N 


3.3(11) 


1.6(11) 


2.1 





5(4) 


HSOj 


3.1(11) 


1.5(11) 


2.1 





3(5) 


S2 


1.4(12) 


2.8(12) 


2.0 





3(5) 


Sect. 4.2: FUV irradiated zone 




Go 


T 










[ISRF] 


[K] 


HSO+ 


3.0(12) 


5.3(11) 


5.6 


106 


101 


ot 


1.3(12) 


2.6(11) 


4.8 


10^ 


550 


CH3OCH3 


2.5(12) 


6.1(11) 


4.1 


10^ 


101 


CH3OH 


5.3(12) 


1.3(12) 


4.0 


lO'^ 


80 


HC3N 


2.9(14) 


7.1(13) 


4.0 


10** 


101 


C2N 


1.0(14) 


2.5(13) 


4.0 


10 


101 



D. SPECIES WITH LARGE DISAGREEMENT 




Go [ISRF] 



Fig. Dl. — Fractional abundance of HSO+ (black solid line) for different FUV fluxes (Go) and FUV attenuation (r). The total density 
was chosen to be lO'' cm~^, while the temperature is fixed at 100 K. The grid points of the curvilinear coordinate system are shown in 
(grey/red) dots and black dotted line. The value of Go (10^) where large deviations were found is indicated by a vertical red/grey dotted 
line. 
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10" 10^ lO"* 10" 

Go [ISRF] 



Fig. D2. — Fractional abundance of O2 (black solid line) for different FUV fluxes (Go) and FUV attenuation (r). The total density 
was chosen to be 10^ cm~^, while the temperature is fixed at 550 K. The grid points of the curvilinear coordinate system are shown in 
(grey/red) dots and black dotted line. The value of Go (10^) where large deviations were found is indicated by a vertical red/grey dotted 
line. 
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